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INTRODUCTION 


Several  Copies  of  importance  for  timely  research  were  undertaken  in 
this  project  this  year.  In  the  first  instance  it  was  found  that 
insufficient  funds  were  available  in  the  Seismic  Research  Center 
Research  contract  (F08606-79-C-0007) ,  Task  4.6)  to  completely  analyze  t* 
for  all  of  the  epicenter-station  combinations  desired.  Thus  funds  from 


the  present  contract  were  used  to  analyze  the  E.  Kazakh-AEDS  paths. 

This  work  was  accomplished  as  discussed  in  the  following  section  and  in 
Appendix  A. 


Another  topic  of  importance  was  the  extension  of  the  data  base  to 
which  the  Generalized  Linear  Model  yield  estimation  work  of  Shumway  and 
Blandford  (1982)  had  been  applied.  Results  from  this  work  are  also 
discussed  in  the  following  section  and  in  Appendix  B. 


There  has  been  considerable  interest  on  the  part  of  workers  in 
Sweden  on  the  possible  uses  of  "dynamic"  criteria  for  improvement  of  the 
automatic  association  program.  For  example,  how  is  one  to  make  use  of 
the  information  that  an  apparently  large  magnitude  event  is  observed  at 
only  a  few  stations  and  is  not  observed  at  some  stations  nearby  the 
hypothesized  epicenter?  Some  work  along  these  lines  is  discussed 
briefly  in  the  following  sections;  and  memos  outlining  the  detailed 
research  are  given  in  Appendix  C. 

Other  related  work  performed  under  this  contract  was  the  generation 
of  synthetic  data  for  testing  of  automatic  association  programs.  The 
techniques  of  generating  this  data  and  the  uses  to  which  it  was  put  are 


discussed  in  an  unsigned  report  written  by  North  and  Olsen  (1983): 

Report  No.  3,  Group  of  Scientific  Experts,  Special  Study  Groups  3  and  5, 
for  the  Conference  of  the  Committee  (CCD)  on  Disarmament.  This  work  is 
also  discussed  in  the  following  section,  and  excerpts  from  the  CCD 
report  are  included  in  Appendix  D. 
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SUMMARY  AND  SUGGESTIONS  poR  FURTHER  RESEARCH 


Kazakh- AEDS  t* 

With  respect  to  determination  of  the  E.  Kazakh  to  AEDS  t*  values 
Appendix  A  gives  the  Abstract,  Table  of  Contents,  Summary  and 
Conclusions,  and  Recommendations  for  Future  Research  from  report 
S!)AC-TR— 82-1  in  which  approximately  150  AF,T)S-E.  Kazakh  t*  values  were 
determined  and  reported.  In  general  terms  the  additional  analysis  did 
not  change  the  qualitative  conclusions  of  the  report;  the  chief  benefit 
of  the  analysis,  of  exclusive  relevance  to  the  AEUS  network,  was  the 
determination  of  the  actual  station  t*  values  themselves  which  can  be  of 
use  in  determining  station  corrections  for  magnitude-yield  and  for 
spectral  discrimination.  The  actual  t*  values  themselves  are  classified 
and  may  be  determined  by  reference  to  the  cited  report. 

General  Linear  Model  Magnitude  Studies 

In  Appendix  R  we  give  a  full  general  linear  model  analysis  of  WWSSN 
SP  data  for  the  events  Shoal,  Piledriver,  Faultless,  Gasbuggy,  Rulison, 
and  Rio  Blanco.  The  analysis  leads  to  the  conclusion  that  a  single 
hard-rock  curve  cannot  encompass  shots  in  granite,  tuff,  and  sandstone- 
shale.  We  find  that  Gasbuggy  lies  0.1  mfc  below  the  Rhoal-P i ledriver 
prediction  for  granite;  and  that  Rulison  and  Rio  Blanco  lies  0.6  m 

D 

below.  It  Is  also  of  interest  that  the  Faultless  explosion  which  has 
been  reported  to  have  a  large  magnitude  for  its  yield  appears  to  have  a 
large  positive  effect  due  to  pP.  When  this  effect  is  removed  the  shot 
lies  near  the  theoretical  curve  passing  through  Shoal  and  Piledriver. 
This  is  of  importance  because  Der  et  al.  (1982)  found  attenuation  under 
Faultless  to  be  comparable  to  that  under  NTS,  and  the  previously 
reported  anomalously  high  magnitude  of  Faultless  has  been  in  conflict 
with  this  observation. 

Dynamic  (Amplitude)  Criteria  for  Automatic  Association 

In  Appendix  C  is  a  copy  of  a  memorandum  outlining  the  work 
performed  to  develop  a  dynamic  method  of  checking  automatic  association 
results.  At  the  time  this  memo  was  written  it  seemed  to  us  that  the 
algorithm  was  ready  to  be  Implemented.  However,  since  that  time  further 
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thought  has  suggested  two  different  possible  modifications  to  the 
procedure . 

For  the  first  possibility,  referring  to  Appendix  C,  in  Figure  1  we 
now  suggest  that,  in  the  uppermost  "box"  that  the  test  at  a  station  be 
whether  what  happened  at  that  station  has  a  small  or  large  probability. 
That  is:  (1)  If  the  signal  is  not  detected,  what  is  the  probability  that 
it  would  be  below  the  noise?  (2)  If  the  signal  is  detected  but  an 
amplitude  is  not  reported,  what  is  the  probability  that  the  signal  is 
above  the  noise?  (3)  If  the  signal  is  detected  and  an  amplitude  is 
measured  what  is  the  minimum  of  the  two  probabilities  (I)  That  it  has  a 
difference  from  the  predicted  amplitude  greater  than  or  equal  to  that 
observed?  (II)  That  it  would  be  detected,  given  its  observed  amplitude? 

The  probability  threshold  would  be  set  low  enough  that  there  would 
be  a  small  chance  of  a  station  or  amplitude  being  thrown  out  in  a  good 
event . 

To  make  these  functional  changes  in  the  subroutines  of  Appendix  C 
would  require  substantial  coding  changes  in  subroutine  detz ,  and 
relatively  minor  changes  in  subroutine  dynkin. 

The  advantages  of  these  changes  are  (1)  The  probability  threshold 
is  more  intuitively  appealing  than  a  scaled  magnitude  threshold.  (2) 

The  program  will  run  much  more  rapidly  because  subroutine  amp  need  only 
be  called  once  in  subroutine  detz  instead  of  ns  (the  number  of  stations) 
times  The  routine  should  run  about  ns  times  more  rapidly. 

The  second  new  overall  approach  would  simply  be  to  modify  the 
existing  "screen"  routine  in  AA  as  discussed  by  Blandford  et  al.  (1983). 
Points  would  be  taken  away  from  the  event  score  if  detections, 
non-detections,  or  observed  amplitudes  were  out  of  line  with  the 
calculated  maximum  likelihood  magnitude.  See  Appendix  C  for  a  second 
memorandum  discussing  this  approach.  At  present  it  seems  to  us  that 
this  approach  is  simpler,  faster,  more  well-integrated  into  normal  AA 
practice;  and  probably  gives  better  results  than  even  the  revised 
approach  of  the  first  memo. 


Further  research  along  these  lines  should  be  to  make  the  above 
■>  re*  r  inn  i  ng  modifications  and  to  apply  and  test  the  techniques  in  an 
operating  AA  program,  perhaps  using  as  data  the  synthetic  data  discussed 
below. 

Synthetic  Data  for  Evaluation  of  Automatic  Association 

In  order  to  provide  a  data  base  for  experiments  in  international 
seismic  communication  using  the  World  Me terological  Network  (WMO)  by 
participants  in  the  Conference  of  the  Committee  on  Disarmament  (CCD); 
and  to  serve  as  a  test  data  base  for  exercising  the  Swedish  and  Center 
for  Seismic  Studies  (CSS)  automatic  association  programs,  it  was  desired 
to  generate  synthetic  arrival  data.  Appendix  D  is  an  extract  from  North 
and  Olson  (1983)  which  describes  fairly  completely  the  technique  of 
generating  the  data.  Here  we  shall  make  a  few  additional  comments. 

The  program  used  to  generate  the  data  was  an  almost  complete 
rewrite  of  a  program  originally  written  by  M.  Chinnery  and  described  in 
Chinnery  (1980).  The  seismicity  distribution  was  retained  unchanged. 
Major  changes  were  made  by  allowing  detection  of  all  phases  to  be 
controlled  by  the  same  detection  criteria  as  the  P  phase  themselves. 

This  required  obtaining  distance-amplitude  curves  for  several  phases  for 
which  such  curves  were  not  available  in  the  literature.  The  curves  were 
ob‘  ined  by  the  use  of  the  results  in  Sweetser  and  Blandford  (1973)  for 
PP  and  all  PKP  branches;  and  analysis  of  phases  measured  by  the  NEP 
analyst  for  the  time  period  26  August  1978,  through  1  May  1979.  The 
phases  analyzed  were  those  with  20  or  more  reported  arrivals  in  this 
time  period,  and  as  seen  in  Appendix  D  included  PcP,  PKKPab ,  PKKPbc, 
PKKPdf,  PKPPKP,  SKPdf,  SKPab ,  ScP,  pP  and  sP.  Scattergrams  were  also 
made  of  the  periods  for  all  phases  as  a  function  of  distance.  The 
phases  SKPab,  PP,  and  PKPPKP  has  periods  of  1.3  +/-  0.5  seconds,  PcP  has 
0.9  +/-  0.5.  All  other  phases  seemed  to  have  periods  of  1.0  +/-  0.5. 

In  the  synthetic  data  these  periods  were  used,  assuming  a  uniform 
distribution  within  the  indicated  limits. 

Point  10  in  Appendix  D  is  incorrect;  the  standard  deviation  of  the 
Gaussian  error  for  each  slowness  component  is  2.0  seconds/degree  for 
small  arrays,  and  1.0  for  LASA  and  NORSAR.  These  are  values  which  we 


have  found  to  be  a  proper  tolerance  in  practice  for  "making"  small 
events.  As  such  they  probably  represent  the  errors  to  be  expected  for 
small  S/N  events;  and  are  much  larger  errors  than  are  found  for  a 
calibrated  array  for  large  events.  Examples  of  such  excellent  location 
capability  may  be  found  in  Dahlman  and  Israelson  (1977).  Typical  errors 
of  200  km  would  correspond  to  a  slowness  error  on  the  order  of  0.1 
sec/km,  one  tenth  the  error  listed  above  for  LASA  and  NORSAR.  In  future 
generation  of  synthetic  data  it  would  be  best  to  make  the  slowness  error 
a  function  of  S/N,  perhaps  varying  from  1.0  to  0.1  as  the  S/N  varys 
from  1.5  to  10. 

As  discussed  in  Appendix  D  there  were  complicated  rules  governing 
the  naming  and  masking  of  detected  phases.  In  retrospect  some  of  the 
conventions  were  not  as  well  thought  through  as  they  should  have  been. 
For  example,  explosions  were  allowed  to  generate  pP  and  sP  phases 
because  insufficient  consideration  was  given  o  the  interference  e f feels 
of  phases  arriving  close  in  time  to  each  other.  A  rule  for  the 
synthetic  data  generation  which  should  have  been  implemented  would  lump 
detected  arrivals  within  3-5  seconds  of  each  other  into  a  single  arrival 
with  the  name  and  arrival  time  of  the  first  phase  and  the  amplitude  of 
the  maximum  amplitude. 

Also,  the  rule  that  phases  were  named  first  and  then  phases  were 
thrown  away  if  masked  by  preceding  phases  occasionally  led  to  anomalous 
results,  related  to  the  fact  that  pP  was  not  named  if  P  was  not 
detected.  Thus,  if  P  was  detected  and  pP  named,  and  if  then  that  P  was 
masked  by  a  preceding  phase  then  the  pP  was  left  in  the  SAQ  without  a  P 
to  go  with,  an  unlikely  analyst  decision.  On  the  other  hand,  perhaps 
this  outcome  could  be  allowed  to  stand  as  an  analyst  error  in  which  the 
analyst  associated  this  phase  as  pP  to  the  masking  phase.  But  further, 
it  is  possible  that  the  masking  phase  itself  was  a  pP.  These  problems 
should  be  worked  out. 

Finally,  we  have  the  suspicion  that  the  later  phase  amplitudes  are 
too  large  in  many  cases.  In  real  life  phases  are  not  detected  if  they 
are  below  the  coda,  and  this  may  badly  bias  the  amplitude  distribution 
as  observed  in  the  NEP  and  other  bulletins.  A  study  should  be  carried 


out  with  the  general  "Ringdal"  approach  to  better  estimate  later  phase 
amplitudes;  always  including  a  coda  measure  for  the  noise  when  a  phase 
is  not  detected.  In  addition,  those  phases  PP,  SKPab,  and  PKPPKP  with 
average  periods  of  1.3  seconds  may  be  harder  to  detect  than  their 
amplitudes  (which  are  corrected  for  period)  would  suggest.  Uncorrected 
for  period,  considering  the  typical  response  of  short-period 
instruments,  their  amplitude  on  the  film  would  be  about  0.5  times  the 
amplitude  of  a  corresponding  1,0  Hz  signal.  Furthermore,  the  noise  at 
that  period  would  be  higher  than  at  1  Hz.  These  problems  should  also  be 
investigated;  and  as  a  quick  fix  we  would  recommend  that  the  amplitudes 
of  these  later  phases  be  reduced  by  1/2  in  any  generation  of  new 
synthetic  data. 
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ABSTRACT  (U) 

(U)  The  values  of  the  apparent  t*  were  estimated  for  a  large  number  of 
paths  from  nuclear  explosions  at  various  test  sites  to  SRO,  AEDS  and  a 
small  set  of  LRSM,  SDCS  and  array  sites.  Multiple  regression  analyses 
of  the  data  showed  that  NTS  explosions  have,  at  teleseismic  stations  a 
t*  .2-. 3  sec  higher  than  similar  explosions  at  the  Kazakh  test  sites. 

The  slightly  higher  (.1  sec)  t*  estimates  for  Yucca  Flats  explosions 
relative  to  Pahute  Mesa  shots  must  be  due  to  differences  in  RDP  possibly 
due  to  the  lower  rigidity  and  residual  gas-filled  porosity  of  the  tuffs 
at  Yucca  Flats,  since  P  wave  spectra  at  near  regional  distances  also 
show  the  same  difference.  The  estimates  of  t*  for  Soviet  PNE  explosions 
in  the  shield  areas  of  the  USSR  are  similar  to  those  derived  for  the 
Kazakh  test  sites.  The  French  nuclear  explosions  in  the  Sahara  show  a 
depletion  in  the  high  frequency  content  in  the  P  wave  spectra  relative 
to  Kazakh  nuclear  explosions  recorded  at  the  same  sites.  The 
contributions  of  the  upper  mantle  structures  under  the  recording 
stations  vary  within  a  few  tenths  of  a  second  in  t*,  WHYK  and  EOBO 
appear  to  be  the  stations  with  the  low  Q  in  the  underlying  mantle. 
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SUMMARY  AND  CONCLUSIONS  (U) 


(U)  The  main  conclusions  of  this  study  regarding  the  variations  of 

attenuation  under  the  source  and  receiver  regions  studied  are  as 

follows ; 

a)  The  most  significant  differences  in  apparent  t*  are  between  the 
source  region  terms  of  Kazakh  and  NTS,  which  are  in  the  range  of  ,2 
-  .3  sec,  with  the  higher  t*  associated  with  NTS, 

b)  There  appears  to  be  a  smaller  but  significant  differential  of  the 
order  of  ,1  sec  in  the  raw  apparent  t*  between  paths  involving  Yucca 
Flat  and  Pahute  Mesa  explosions  at  common  stations.  A  possible 
source  function  with  low  overshoot  (B  -  0)  and  lower  rise  time  due 
to  the  soft  sediments  at  Yucca  Flats  may  be  the  explanation. 

c)  The  contributions  of  the  upper  mantle  under  the  receiving  stations 
show  smaller  relative  variation  than  the  test  site  (source  region) 
terms.  Stations  with  relatively  higher  attenuation  in  the  upper 
mantle  are  WHYK  and  ZOBO,  The  relative  variations  of  the^station 
terms  remain  similar  after  the  source  spectra  are  factored  out. 

d)  The  available  data  from  the  French  nuclear  explosions  in  the  Sahara 
indicate  a  significant  loss  of  high  frequency  content  in  the  P  waves 
relative  to  those  originating  at  Kazakh.  This  is  indicative  of  high 
attenuation  under  the  Sahara  test  sites, 

e)  The  t*  estimates  from  Soviet  PNE  explosions  in  shield  areas  are  not 
significantly  different  from  those  obtained  from  nuclear  explosions 
at  the  two  main  Kazakh  test  sites  at  common  stations, 

f)  The  spectra  of  the  two  WUS  granite  explosions  PILEDRIVER  and  SHOAL, 
when  scaled  to  a  common  yield  at  the  common  station  NPNT ,  are  not 
significantly  different  from  those  of  other  NTS  explosions,  but  show 
less  high  frequency  content  than  the  Kazakh  explosions, 

g)  Comparison  of  spectra  for  common  events  recorded  at  RKON  and  NORSAR 
indicate  a  slightly  higher  attenuation  in  the  mantle  under  RKON 
compared  to  NORSAR,  More  data  is  needed,  however,  to  establish  this 
with  certainty, 

h)  Analysis  of  short  period  S  waves  from  shallow  earthquakes  in  plates 
or  at  the  edges  of  the  Russian  shield  that,  even  assuming  a  delta 
function  source,  the  waveforms  are  inconsistent  with  an  apparent  t* 
of  4  sec;  the  best  fit  is  around  2  sec.  The  values  of  the  observed 
S/P  ratios  may,  on  the  other  hand,  be  reconciled  by  higher  t* 
(absolute) , 

Other  results  of  this  study  are  summarised  below: 


a)  Most  P  waves  from  the  nuclear  explosions  analysed  contained  source 
related  secondary  arrivals  besides  the  surface  reflection  pP, 

b)  Spectral  ratios  of  pairs  of  events  at  many  sensors  of  NORSAR  become 
dissimilar  beyond  2-3  Hz,  indicating  that  the  concept  of  "relative 
receiver  functions"  breaks  down  at  high  frequencies, 

c)  Determination  of  the  times  and  amplitudes  of  secondary  arrivals  is 
inherently  limited  in  accuracy  by  the  uncertainties  of  the  source 
time  function  (self-noise)  and  useful,  coherent  signal  bandwidth, 
which  is  in  turn  determined  by  t*,  background  noise  and  random 
fluctuations  in  the  media  of  propgation.  These  factors  can  be 
evaluated  for  each  event,  and  the  uncertainties  in  the  estimated 
source  depths  and  yields  can  be  formulated  in  terms  of  bandwidth  and 
signal  to  noise  ratio. 


RECOMMENDATIONS  FOR  FUTURE  RESEARCH  (U) 


(U)  Although  the  results  presented  in  this  report  give  a  good  first 
look  at  the  variations  of  Q  under  some  important  test  sites  and 
recording  stations,  it  is  clear  that  the  amount  of  data  used  was  not 
sufficient  to  determine  recording  station  t*  terms  with  the  desired 
accuracy  of  ,05  sec  in  95T  confidence  Limits,  Moreover,  some  maior  data 
sets,  such  as  the  AEDS ,  SDCS  and  SRO  data,  do  not  have  the  necessary 
overlap  to  determine  the  relative  t*  values  among  stations  in  the 
various  data  sets  accurately.  The  coverage  of  the  various  test  sites  by 
the  SRO  stations  is  also  non-overlapping  and  sparse.  There  is  a  need, 
therefore,  to  expand  the  data  sets  and  provide  more  overlap  with  common 
events  at  common  test  sites  as  far  as  the  limitations  of  geography 
allow.  The  AWRE  array  data  can  be  extremely  useful  in  providing  both 
the  overlap  and  standards  of  common  reference  because  of  the  Long 
operational  duration  of  these  arrays  which  overlap  in  time  all  of  the 
other  data  sets.  Early  in  this  study,  we  decided  that  the  hand- 
digitized  WWSSN  data  is  not  likely  to  be  reliable  for  measurement  of  t*, 
which  critically  depends  on  high  frequency  information  (as  we  have 
emphasized  repeatedly  in  this  report).  We  plan  to  study  these  stations 
by  using  short  period  S  wave  data  from  deep  earthquakes  which  were 
proven  to  be  sensitive  indicators  of  t*  variation  in  past  studies  (Der, 
McElfresh  and  O'Donnell,  1982;  Der,  Smart  and  Chaplin,  1980),  In  doing 
this,  we  plan  to  make  extensive  use  of  the  WWSSN  chips  available  at  the 
Center  for  Seismic  Studies. 

(U)  More  work  needs  to  be  done  to  explore  the  fundamental  limitations 
on  the  information  about  the  source  and  path  properties  contained  in 
band-limited  data.  In  the  past,  the  accuracy  of  the  results  obtained 
was  considerably  over  estimated,  or  the  information  available  was  not 
optimally  utilized.  Clearly,  it  is  futile  to  attempt  to  derive  source 
information  beyond  some  fundamental  limits  imposed  by  the  data. 
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APPENDIX  B 


Blandford  and  Shuraway  (1982),  hereafter  referred  to  as  paper  A, 
developed  a  Generalized  Linear  Model  (GLM)  technique  for  simultaneously 
estimating  event  magnitudes  and  station  corrections  while  also  retaining 
the  advantages  of  the  Ringdal  (1976)  maximum  likelihood  magnitude 
estimation  procedure  in  which  bias  (due  to  the  fact  that  small 
amplitudes  are  not  detected  and  thus  not  reported)  is  minimized. 

The  major  task  undertaken  in  the  present  study  on  the  subject  of 
magnitude :yield  is  the  analysis,  using  the  GLM,  of  an  additional  4 
off-NTS  underground  explosions  in  the  Southwest  United  States.  These 
events  offer  us  excellent  opportunities  for  determining  true  errors  in 
magnitude  .’yield  estimation;  and,  as  we  shall  see,  the  application  of  the 
GLM  is  important  for  3  out  of  the  4  events  because  their  true  magnitudes 
are  less  than  5.0. 

In  Table  la  we  see  the  6  events  considered  in  the  present  study, 
Rulison,  Gasbuggy,  Faultless,  and  Rio  Blanco,  together  with  two  events 
previously  analysed  in  A:  Piledriver  and  Shoal, 

The  events  were  analyzed  in  this  study  using  exactly  the  same 
procedures  as  in  A.  Useful  readings  were  obtained  from  95  WWSSN 
stations,  and  the  number  of  signal  detections,  as  distinct  from  noise 
measurements,  ranged  from  71  for  Faultless  to  13  for  Gasbuggy. 

In  Table  lb  we  see  the  results  of  the  maximum  likelihood 
calculations.  The  maximum  likelihood  magnitudes  for  the  smaller  events 
are  0.34  to  0.60  m^  smaller  than  the  simple  average  of  the  detecting 
stations,  and  also  are  smaller  than  the  simple  magnitudes  calculated  by 
Marshall  et  al.  (1978).  This  latter  comparison  suggests  that  there 
were  signals  below  the  noise  in  the  LRSM  network  from  which  detection 
were  selected  by  Marshall  et  al. 

Note  that  the  Faultless  maximum  likelihood  magnitude  is  6.49  as 
compared  to  the  simple  average  of  6.40;  this  reflects  the  presence  of  9 
clipped  signals  out  of  80  detections. 


Date  Time  Location  Depth(m)  Yield  (kt)  Medii 


Rulison  I  69/09/10  21:00 


Gasbuggy  67/12/10  19:30 


Fault-  I  68/01/19  18:15 


Blanco 


73/05/17  16:00 


e-  66/06/02  15:30 

Driver 


Shoal  I  63/10/26  17:00 
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36.68 
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3.26  +.00 

4.7'> 

Gasbuggy 

5.07 

29 
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3.14 

3.13  -.01 
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( e  s 

600 
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.80  (4) 
(.90)  (5) 

4.14 

4.52  + . 38 

6.33 

Rio 

Blanco 

4.80 

5.17 

5.31 

90 

1.30  (3) 
(1.15)  (6) 

3.56 

3.56  +.00 

5.02 

Pile 

Driver 

5.49 

5.57 

5.71 

62 

.25,  (.17) 

3.42 

3.64  +.22 

5.49 

Shoal 

4.75 

5.09 

5.01 

12 

•  21,  (.15) 

2.7  7 

2.99  +.22 

4  .  /  5 

(1)  Marshall  (1972) 

(2)  Cohen  (1970) 

(3)  Marshall  et  al .  (1978) 

(4)  Springer  (1974) 

(5)  Frazier  (1972) 

(6)  von  Seggern  (1974) 

(7)  Lg  ratio  to  Greeley  at  common  good  quality  stations,  See  Clark  (  19bH> 


m^  -  maximum  likelihood 

Tn^  -  simple  average  of  detecting  stations 
mMSR  _  from  Marshall  et  al.,  (1978) 
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Theoretical  seismograms,  shown  in  Figure  1,  were  calculated  for 
each  of  the  events  using  the  yield  and  pP  delay  parameters  given  i i 
Table  lb,  a  pP  reflection  coefficient  varying  according  to  0.5(l+exp 
(~f**2))  and  a  t*  value  of  0.4,  taken  to  be  typical  of  the  average  path 
to  a  detecting  station.  The  difference  between  log(A/GT),  where  G  is 
the  instrument  correction  at  period  T,  for  the  waveform  without  pP  and 
the  waveform  with  pP  is  defined  as  the  pP  "correction"  and  is  subtracted 
from  the  observed  m^  value  for  the  corresponding  event  in  order  to  take 
out  the  pP  effect. 

As  in  A,  however,  the  Piledriver  pP  correction  is  added  back  in  for 
all  events.  In  this  way  the  actual  Piledriver  m,  calculated  is  retained 

D 

as  a  standard  and  all  other  events  are  adjusted  relatively  to  it.  Thus 
there  results  an  m^:yield  curve  where  the  effects  of  pP  have  been 
removed,  but  the  Pildriver  m^  has  the  same  numerical  value  as  before  the 
pP  correction. 

In  calculating  the  pP  corrections  the  pP  times  used  were  the 
preferred  one  which  are  not  enclosed  in  parentheses  in  Table  I. 

Figure  2  shows  the  results  of  these  calculations.  For  those  shots 
in  granite  and  tuff,  the  magnitudes  corrected  for  pP  are  seen  to  fall 
close  to  lines  with  the  predicted  average  theoretical  slopes;  1.0 
between  Shoal  and  Piledriver,  0.7  between  Piledriver  and  Faultless. 

However,  the  other  shots,  in  shale  and  sandstone  are  seen  to  fall 

well  below  this  line,  (although  the  two  events  within  56  kilometers  of 

each  other,  Rulison  and  Rio  Blanco  fall  within  0.05  m,  of  a  line  with 

b 

slope  1.0)  suggesting  that  a  curve  which  is  adequate  for  granite  and 
tuff  might  not  be  adequate  for  sandstone  and  shale.  This  suggests  that 
a  single  "hard-rock"  curve  could  not  be  suitable  for  shots  in  such 
different  media. 

Note,  however,  that  if  the  maximum  likelihood  estimates  for 
magnitude  has  not  been  made,  all  the  events  would  have  fit  fairly  well 
on  a  single  curve  of  lesser  slope. 


log  Y 


The  observed  discrepancy  between  the  Piledriver  and  Oasbuggy  m^ 

values  is  consistent  with  observed  reduced  displacement  potentials 

(RDPs)  for  this  event.  Examination  of  Figure  la,  b;  from  Murphy  (1978) 

and  Murphy  and  Bennett  (1980)  suggests  that  the  Oasbuggy  POP  would  be 

8 

between  6000  and  8000  m  ,  while  the  Piledriver  ROP  would  be  between 

20000  and  80000.  Thus  we  might  expect  an  m  difference  around 

b 

60000/6000  or  1.0  m^  units,  with  a  range  from  70000/6000  to  80000/6000 

or  0.52  to  1.17  m  .  The  observed  difference  is  0.61  m  ;  consistent  with 

b  b 

the  ROP  differences  and  much  greater  than  the  0.7  m^  expected  from  the 
Shoa 1-P i ledr iver  curve. 

Consistency  with  the  observed  ROP  values  suggests  that  the  even 

lower  scaled  m,  values  for  Rulison  and  Rio  Blanco  are  not  a  mistake  but 
b 

reflect  true  differences.  Note  that  even  these  0.6  m,  low  values  are 

b 

within  the  range  allowable  by  the  Piledriver/Oasbuggy  ROP ' s . 

Tn  this  respect  it  is  puzzling  that  Nuttli  (personal  communication, 
1987)  found  excellent  agreement  between  T.  magnitudes  and  yields  for 

n 

these  events.  In  Figure  4  we  see  that  there  is  very  small  scatter  in 

the  magnitude  yield  relationship.  This  is  consistent  with  earlier  work 

by  Blandford  and  Klouda  (1980)  and  others  showing  very  weak  dependence 

of  magnitude  on  medium;  but  seems  to  be  inconsistent  with  the  large  ROP 

differences.  Perhaps  the  ROP ' s  are  not  reliable;  we  now  know  that  they 

are  measured  in  the  non-linear  regime,  and  are  thus  not  truly 

representative  of  the  far-field.  Thus  perhaps  here  again  we  see  the 

impressive  capability  of  from  to  give  a  magnitude  independent  of 

medium  while  m,  from  P  waves  shows  a  large  effect, 
b 

We  see  in  Figure  2  that  after  correction  for  pP  Faultless  is  close 
to  the  theoretical  curve.  The  pP  correction  for  Faultless  is  very 
large,  0.78  m.  units  as  compared  to  0.22  for  Piledriver.  Most  of  this 

D 

correction  can  be  traced  to  the  unusually  long  period  of  the  synthetic 
Faultless  seismogram,  1.7  seconds,  which  can  itself  be  traced  to  the 
fortuitous  timing  of  the  pP  signal.  This  can  be  seen  by  comparison  of 
Piledriver  and  Faultless  in  Figure  1.  That  this  effect  is  operative  in 
the  actual  observations  was  checked  by  comparing  periods  at  common 
stations  for  these  two  events.  The  values  for  Piledriver  and  Faultless 


were  0.9  and  1.4  respectively  as  lompared  to  the  theoretical  values  of 
0.7  and  1.7  as  seen  in  Figure  1.  (Note  that  the  theoretical  difference 
without  pP  is  only  0.8  and  1.0).  The  fact  that  the  difference  is  in  the 
proper  direction  and  general  magnitude  is  encouraging.  The  discrepancy 
in  precise  periods  could  easily  be  due  to  analyst  bias  toward  periods  of 
1.0  and  should  be  further  investigated  by  careful  analysis  at  stations 
which  possess  excellent  signal  definition  and  at  digital  stations. 

The  fact  that  after  correction  for  pP  Faultless  falls  close  to  a 
theoretical  curve  suggests  that  signals  emerging  from  the  Faultless  test 
sites  are  not  anomalous  as  compared  to  those  emerging  from  NTS.  Thus 
the  fact  that  Der  et  al.,  (1982)  found  t*  at  the  Faultless  site  to  be 
close  to  that  at  NTS  does  not  seem  to  be  anomalous.  On  the  other  hand, 
it  seems  that  the  fact  that  amplitudes  of  received  signals  at  the 
Faultless  site  average  0.3  m^  low  (Der  et  al.,  1982)  is  still  a  paradox. 
Reciprocity  would  suggest  that  the  Faultless  m^  should  fall  below  the 
curve  in  Figure  2a;  not  right  on  or  slightly  above  as  it  does. 

As  a  secondary  topic  in  this  study  we  mention  that  in  paper  A  a 
literature  survey  was  made  to  determine  the  geological  layer  parameters 
of  the  WWSSN  stations  used  in  the  analysis;  and  Haskell-type 
calculations  were  made  to  determine  the  resulting  elastic  station 
correc t ions . 

As  one  task  in  the  present  work  we  have  tested  whether  application 
of  these  station  corrections  appears  to  offer  any  benefit;  and  the 
conclusion  is  that  apparently  it  does  not.  (The  conclusion  is  not  too 
surprising  the  rms  station  correction  is  0.05  m^  units;  while  the  rms 
residual  is  approximately  0.3  m^  units).  In  particular  the  maximum 
change  in  an  ra^  estimate  i6  0.01,  and  the  rms  residual  changes  from, 
0.312  to  0.314;  actually  an  increase! 
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MEMORANDUM 


A.  Kerr,  W.  Dean,  J.  Goncz,  R,  North,  R.  Slunga,  R.  Shumway 


FROM: 


R.  Blandford 


SUBJECT:  Dynamic-Kinematic  Criteria  for  Event  Reality 


DATE: 


7  October  1982 


We  now  have  a  complete  "system"  for  dynamic-kinematic  amplitude 
checking;  attached  to  this  memo.  I  think  this  system  could  be  checked 
out  fairly  well  with  the  artificial  data  John  Goncz  and  I  are 
generating.  The  subroutine  amp  may  diverge  if  very  large  or  small  ID  * 

1  type  values  are  in  the  data.  Bob  Shumway  and  I  will  soon  distribute 
an  improved  version  of  amp,  I  think  this  code  is  in  a  "final"  state 
suitable  for  use  in  existing  AA  programs.  Only  the  I/O  needs  to  be 
modified,  and  this  is  fairly  well  isolated  in  the  event. f  and  mkpkin.f 
files. 

As  long  as  we  are  talking  about  AA,  I  propose  a  measure  of  AA 
quality,  "Association  Percentage"  defined  as  follows.  Define  events 
which  are  to  be  detected;  e.g.  events  with  5  non-array  stations;  plus 
events  with  two  arrays  and  one  confirming  station,  plus  events  with  2  P 
waves  and  one  S,  etc.  All  arrivals  in  this  set  of  events  when  perfectly 
analyzed  constitute  the  "base  set  of  arrivals".  Now  we  go  down  the 
perfect  bulletin  and  pick  out  the  first  event.  Is  there  an  event  in  the 
AA  with  (75%)  of  the  arrivals.  (Obviously  75%  is  a  variable;  probably 
anywhere  between  60%  and  90%  is  alright.)  If  so,  then  all  correct 
arrivals  in  that  AA  event  may  be  added  to  the  total  correct  arrivals 
count;  and  all  incorrectly  associated  arrivals  should  be  subtracted.  I 
propose  that  a  correctly  associated  arrival  should  count  even  if  the 
phase  name  is  incorrect.  (Or  perhaps  it  should  get  1/2  or  0  weight  if 
incorrectly  identified.)  In  this  way  all  "perfect  bulletin"  events  are 
considered.  Any  remaining  arrivals  which  have  been  associated  into 
events  (which  are,  presumably  "false"  or  "bad"  events  and  which  will 
give  an  analyst  severe  trouble  to  clean  up)  will  be  subtracted  from  the 
total  correct  arrivals  count.  The  Association  Percentage  is  calculated 
as  (total  correct  arrivals  count)/(base  set  of  arrivals). 

This  measure: 

Doubly  penalizes  the  creation  of  false  events  and  incorrect 
association  to  good  events. 

Rewards  correct  identification  of  later  phases. 

Does  not  depend  on  location  error  which  is  highly  variable  depending 
on  event  size  and  travel-time  residuals. 

Can  be  easily  calculated  automatically. 

Perhaps  too  severely  penalizes  splitting  large  events,  especially  an 
equal  split. 
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The  problem  of  verifying  that  magnitude  estimators  for  events  which 
have  been  created  by  an  automatic  association  (AA)  routine  are  base  1  u 
a  reasonable  subset  of  detecting  and  non-detecting  stations  has  been 
considered  by  ELvers  (1980).  She  developed  a  "plausibility  function" 
for  the  amplitude  distributions  which  is  essentially  proportional  to  the 
likelihood  of  observing  a  particular  configuration  of  amplitudes  when 
the  event  is  assumed  to  have  occurred. 

In  this  discussion,  we  attempt  to  improve  several  aspects  of  this 
procedure.  As  a  computational  improvement  we  also  have  developed  a 
Newton-Raphson  procedure  for  computing  the  maximum  likelihood  magnitude, 
assuming  the  signal  and  noise  variances  are  known.  This  procedure  runs 
much  faster  than  the  search  techniques  used  by  previous  authors. 

We  also  give  a  method  for  discarding  single  stations  with  outlying 
amplitudes  (I,  see  Figure  1)  based  on  the  "influence"  they  exert  on  the 
magnitude  estimator.  This  replaces  the  Elvers  procedure  of  comparing 
the  separate  components  of  the  likelihood  for  the  purpose  of  determining 
the  influential  stations.  Our  procedure  for  discarding  stations  is  also 
apparently  not  so  likely  to  "destroy"  the  event  in  the  case  that  a  large 
number  of  stations  are  "down"  and  fail  to  report,  but  do  not  inform  the 
bulletin  analysis  center  of  this  fact. 

Slunga  (personal  communication,  1982)  has  pointed  out  that  the 
"dynamic"  amplitude  criteria  of  Elvers  seem  to  erroneously  discard 
arrivals  from  an  event  on  the  basis  of  amplitude  anomalies.  On 
occasion,  although  the  amplitudes  are  indeed  anomalous,  the  probability 
that  the  arrival  time  could  agree  by  accident  seems  even  smaller,  so 
that  it  seems  unreasonable  to  totally  discard  the  arrival  from  the 
event . 

To  help  remedy  this  paradox  we  propose  to  measure  the  overall  event 
plausibility  with  a  "kinematic"  likelihood  ratio  (II)  in  parallel  with 
the  dynamic  one.  If  the  dynamic  criteria  suggest  that  an  observed 
amplitude  is  unreliable  but  the  kinematic  criteria  are  favorable,  then 
we  suggest  that  the  appropriate  procedure  in  most  cases  is  to  assume 


chat  there  is  some  error  in  the  amplitude  measurement  and  to  recalculat 
tho  in.tgn  i  t  ml**  ,  after  making  an  appropriate  change  in  t  li  t*  amplitude  data 
but  to  continue  to  use  the  reported  arrival  time  for  estimation  of 
location. 

Finally,  we  have  developed  a  procedure  for  tying  these  new 
computational  techniques  together.  In  general  terms  this  procedure  is 
not  derived  from  theoretical  considerations  but  embodies  common  sense 
ideas  about  the  causes  of  erroneous  amplitude  measurements. 

It  is  important  to  emphasize  that  these  techniques  cannot  simply  be 
inserted  into  an  automatic  association  program  as  mathematical  routines 
like,  for  example,  the  cosine  function.  Instead  the  routines  contain 
parameters  which  are  specific  to  the  network,  detectors,  and  methods  of 
analysis.  These  parameters  must  be  determined  by  careful  statistical 
analysis  of  network  bulletins  which  are  relatively  error-free,  and  by 
analysis  of  false  events  produced  by  the  automatic  association  program 
in  which  they  are  to  be  used. 


FUNCTIONAL  DESCRIPTION 


In  Figure  i  we  see  a  flow  chart  of  the  proposed  proceduie.  1  u._ 

Roman  numerals  I-II  denote  the  computational  techniques  referred  to  in 
the  Introduction. 

At  the  top  of  Figure  1  we  see  that  the  event  is  submitted  to  the 
program  for  analysis.  Obviously  only  those  events  should  be  accepted 
for  which  the  program  is  prepared.  For  example  we  may  anticipate  that 
the  first  version  of  the  program  will  be  prepared  to  work  only  on 
initial  P  waves.  If  the  event  is  made  up  of  Pg ,  S,  Lg ,  etc.  then  the 
algorithm  has  nothing  to  contribute.  On  the  other  hand  an  event  with 
only  one  P  wave  and  without  an  amplitude  reading  may  be  submitted.  The 
single  detection  might  be  unlikely  because  many  stations  with  lower 
thresholds  nearby  the  detecting  station  did  not  detect.  When  this  P 
wave  is  discarded  then  the  event  might  be  unreliable  on  kinematic 
grounds  . 

The  event  data  is  then  submitted  to  a  subroutine  (I)  which  computes 
the  Ringdal  magnitude  and  detects  outlying  magnitudes.  *  In  Appendix  I  is 
a  discussion  and  liscing  of  a  subroutine  (amp)  which  rapidly  computes 
the  Ringdal  magnitude  given  station  amplitudes  and  thresholds  in  terms 
of  magnitudes,  and  assuming  that  the  signal  and  noise  variances  are 
known.  Before  this  subroutine  can  be  used,  of  course,  distance- 
amplitude  relations  must  be  used  to  transform  amplitudes  and  thresholds 
to  magnitudes. 

The  amplitude  patterns  can  be  subdivided  roughly  into  n^  observed 
amplitudes,  n^  above  the  noise  threshold  and  with  arrival  times  but  not 
amplitudes  reported,  and  n^  below  the  noise  threshold  and  not  detected. 
Of  course  other  subcategories  such  as  clipping  might  be  added.  For  the 
general  unknown  magnitude  case  the  log- 1  ike  1 ihood  as  given  by  Elvers 
(1980)  would  be  of  the  form: 


FLOW  CHART  FOR  DVNAMIC-KINEMATIC  ANALYSTS  FOR  AUTOMATIC  ASSOCIATION 


Enter  Event 


"Large"  effect  on 
m(,  for  some  station?  — 


(PkH  real)  <0.95 
II 


Determine  Event  Type 


Detecting  Station? 


Discard  Station  -  Amplitude  Measured? 


Change  Station 
type  to  ">  noise' 


Pa  -  Amplitude  probability 

(G)  -  Good  event  exit 

Pk  1  -  Kinematic  probability  with  one  less  station 
if  a  suspect  station  has  been  removed 

(B)  -  Bad  event  exit 

Figure  1 


where  m^  are  the  observed  amplitudes,  m  is  the  theoretical  magnitude 


and  : 


m  -  D .  -  C 
J 


ZJ  = 


/ 


+  a 


j 


is  a  standard  value  depending  on  the  mean  noise  level  D.  with  variance 


2  2 
q j  ,  the  theoretical  magnitude  m  with  observed  variance  og  ,  and  a 

signal-to-noise  constant  C.  Since  m  is  the  only  unknown,  (1)  can  be 


maximized  using  the  Newton-Raphson  or  scoring  algorithms,  which  leads  to 


the  maximum  likelihood  estimator  m  and  its  estimated  variance  o~  (see 

m 

Appendix  I).  The  above  approach  is  essentially  that  followed  by  Elvers 


(1980)  . 


As  a  modification  to  the  above  Elvers  (1980)  suggests  comparing  the 
value  of  the  maximized  log  likelihood,  say  In  L(m)  with  some  arbitrary 
threshold,  with  the  first  component  of  (1)  modified  to  behave  like  an 
estimated  interval  probability.  If  the  total  plausibility,  say  In  L(m) 
is  too  small,  then  the  event  can  be  rejected.  The  other  decision  that 
can  be  made  at  this  point  in  the  Elvers  procedure  is  to  discard  a  single 
station  value  based  on  a  single  station  component  of  In  L(in)  .  For 
example  the  decision  to  reject  a  station  j  which  does  not  observe  is 
based  on  the  value  of  ln($(-Zj)),  which  is  the  probability  of  observing 
a  value  below  the  noise  threshold  at  the  jth  station. 

The  problem,  mentioned  by  Elvers  (1980),  with  this  approach  is  that 
the  total  plausibility  depends  in  part  upon  point  values  of  the  density 
and  in  part  upon  interval  probabilities  which  are  integrated  densities. 
Thus  it  is  likely  that  eliminating  an  event  purely  on  this  criterion 
will  unnecessar i ly  eliminate  events.  The  problem  is  magnified  when  one 
uses  the  Elvers  (1980)  procedure  to  eliminate  single  stations;  purely  on 
the  basis  of  their  plausibility  values. 

As  a  better  method  for  detecting  station  outliers,  we  suggest 
recalculating  the  magnitude  estimator  with  each  station  value  missing  to 
obtain,  say  m^_^  for  i“l,...,n,  where  m^  ^  denotes  the  magnitude 
estimated  when  station  i  is  left  out  of  the  calculation.  These  can  be 
compared  with  the  original  estimator,  say  ra,  to  determine  the  "influence 


of  station  i.  For  example,  a  common  measure  of  influence  is  "Cook's 
distance",  defined  in  this  case  by: 


An  heuristic  procedure  suggested  by  Cook  is  to  reject  station  i  when 
z^  .  ^  moves  further  away  from  zero  than  z(a/2),  where  a  is  some 
arbitrary  probability  value. 

As  an  example  of  how  the  two  procedures  for  amplitude  verification 

might  work  in  practice,  consider  the  11  station  "event",  shown 

schematically  in  Table  1.  We  note  first  that  taking  the  simple  mean  of 

the  observed  magnitudes  yields  "m  ■  4.04,  whereas  the  maximum  likelihood 

estimator  is  m  *  3.78  with  estimated  standard  deviation  0"  =  .12  so 

m 

that  the  censored  observations  definitely  pull  the  estimator  down. 

In  order  to  determine  whether  any  of  the  station  values  m^  are 
outliers,  we  may  use  the  analogue  of  Cook's  distance  given  in  equation 
(2).  The  estimated  deleted  station  magnitudes  for  i*l,2,...,ll 

are  shown  in  Table  1  and  we  note  that  the  greatest  change  is  produced  by 
omitting  station  6(m_g  *  3.89)  which  is  the  value  known  to  be  below  its 
noise  threshold  "D^  *  3.  The  elimination  of  this  station  then  moves  the 
maximum  likelihood  estimator  back  towards  the  mean. 

A  program,  DETZ,  which  produces  the  values  of  is  given  in 

Appendix  I. 


TREATMENT  OF  STATIONS  WITH  ANOMALOUS  AMPLITUDES 

Returning  to  Figure  1  we  must  ask  what  course  of  action  to  take  when 
an  anomalous  station  amplitude  is  discovered  by  procedure  I.  First  we 
determine  if  the  most  anomalous  station  has  reported  a  detection  or  if, 
as  in  Table  1,  it  has  simply  failed  to  report  at  all,  implying  that  the 
signal  level  is  below  the  noise  level.  In  the  latter  case  we  follow  the 
approach  suggested  by  Elvers  (1980)  and  simply  discard  the  station. 
However,  we  may  look  ahead  a  bit  and  note  that,  unlike  Elvers  we 


TABLE  1 


Case  I 


Scat  ion 


Data 


m  . 

o. 

InL 

-i 

m 

3.74 

.13 

1.03 

3.80 

.12 

.96 

3.67 

.13 

2.27 

3.74 

.13 

1.03 

3.70 

.13 

1.51 

3.89 

.17 

6.14 

3.80 

.12 

3.06 

3.79 

.12 

2.87 

3.81 

.12 

3.20 

3.78 

.12 

2.73 

3.78 

.12 

2.68 

.  32 
-.25 
.89 
.32 
.60 


-.  14 
-.29 
-.05 
-.00 


immediately  test  for  the  kinematic  probability  of  this  event  (II)  and 
if  Lfio  probability  of  the  event  being  rear  is  less  than  say  0,9.  ih.i  w. 
exit  the  routine  with  a  bad  event  flag.  That  is,  non-detection  by  this 
station  suggests  that  there  is  something  wrong  with  the  event,  and 
unless  the  kinematic  probability  is  high  the  event  is  discarded.  On  the 
other  hand,  if  the  kinematic  probability  of  the  event  is  high,  then 
probably  there  was  something  wrong  at  the  station,  and  we  return  to 
further  analyze  the  event  without  considering  this  station.  Since  in 
practice  many  stations  fail  to  report  due  to  operational  consideration 
we  should  probably  set  a  less  restictive  threshold  for  the  non-detecting 
stations,  that  is,  we  should  tend  to  reject  them  first.  For  example,  if 
for  a  non-detecting  station  z^  ^  -0.7  reject  it  (see  Table  I)  but 

otherwise  require  z^_^  1,0  for  rejection.  Probably  we  should  have 

no  anomalous  non-detecting  station  remaining  before  considering  other 
stations.  Further  work,  perhaps  empirical  is  needed  to  best  determine 
those  thresholds. 

If  the  anomalous  station  did  detect  we  ask  if  an  amplitude  was 
measured  or  not.  If  no  amplitude  was  measured  then  we  must  have  the 
situation  where  this  station  could  not  have  been  expected  to  detect; 
typically  a  station  in  a  shadow  zone  to  a  small  event.  Again  we  discard 
the  station  and  check  the  kinematic  probability.  Since  we  have  by  this 
time  eliminated  many  of  the  incorrectly  non-reporting  stations,  the 
event  magnitude  should  not  be  biased  too  low  so  that  stations  should  not 
be  incorrectly  thrown  out  at  this  point. 

Finally  we  have  the  case  where  there  _i^  an  amplitude  measurement  at 
the  anomalous  station:  either  it  is  too  large  or  too  small.  In  either 
case  we  assume  that  some  blunder  has  been  made  and  that  if  a  true  signal 
has  been  detected  then  the  amplitude  has  been  recorded  or  transmitted 
incorrectly.  Thus  we  simply  change  the  amplitude  measurement  (not  in 
the  original  files  of  course)  to  state  that  the  signal  was  greater  than 
the  noise.  Again  we  proceed  through  the  kinematic  criterion.  Note  that 
if  this  event  comes  around  again  then  this  detection  will  not  have  a 
measured  amplitude  so  that  if  it  is  unlikely  that  this  station  could 
have  detected  the  event  at  all  the  arrival  will  be  discarded. 
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KINEMATIC  REJECTION  CRITERIA  (II) 


The  kinematic  criteria  for  the  event  e  *  i  •,  t  n .  .  u  i 
as  follows.  First,  it  is  necessary  to  establish  event  types  which  will 
consist  of  the  number  of  array  stations  and  the  number  of  non-array 
stations  detecting.  These  stations  may  also  have  to  be  broken  down  by 
the  criteria  of  analyst  and  automatic  detection.  For  each  type  of  event 
we  may  determine  the  ratio  of  good  events  to  false  events  of  each  type 
in  each  day.  This  ratio  then  gives  directly  the  probability  that  an 
event  of  each  type  is  real. 


PROCESSING  AFTER  DYNAMIC  CHECKING 


If  any  detections  have  been  discarded  then  the  reduced  set  of 
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3  In  L(m) 


n2  (z  R(Z,)+:r(Z,))  "3  (Z.R(Z. )+R  (-Z. ) ) 


Then  the  Newton-Raphson  interations  are  of  the  form 


/  3  In  L(i.i)  v 
dm  mi 


(t)  In  L(m)  ^ 


with  the  final  estimator  m  having  an  estimated  standard  error  given  by 


S2ln  L(m)v% 

"  *  <*,*  ). 

m 
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c 

c...  Driver  program  for  dynamic-kinematic  subroutine 

c 

Integer  DIMEN, maa, mad, msa,msd 
Parameter  (DIMEN*50) 

Parameter  (maa*5,mad*0,msa*7,msd*0) 

Real  ym(DIMEN) ,no(DIMEN) ,sg(DIMEN) 

Real  ymc (DIMEN) 

Real  f  m, se , tl , as , c , znop ,  zlrg , zv lrg , pthr , pobs 
Real  pkin(0:maa,0:madt0:fflsai0:msd) 

Integer  id (DIMEN) 

Integer  idc(DIMEN) ,gb (DIMEN) ,gbc (DIMEN) , ns, i,naa,nad,nsa,nsd 
Integer  j,k,l 

Logical  ary (DIMEN ) ,anlyst (DIMEN) 

Logical  ok 
c 

c...  See  comments  in  dynkin  for  the  following  parameters 

c 

88*. 4 
c*l . 

znop*-. 7 
zlrg*+l . 
zvlrg*l .5 
pthr*. 90 

naa*maa 

nad*mad 

nsa-msa 

nsd*msd 

c 

c...  Get  the  kinematic  matrix 

c 

call  nkpkin(pkin,naa,nad,nsa,nsd) 

open  (2, statu s*'new'  , form* ''formatted' ,file*'event. out') 
c  write(2,p0)((((pkin(i, j ,k,l),i, j,k,l,i-0,naa), j*0,nad),k*0,nsa 

c  &  l*0,nsd) 

50  format(fl0.2,4i5) 
c 

c...  Get  event 

c 

call  event(ym,no,ns,id,sg,ary,anlyst,gb) 
c 

vrite(2,l00)  (i,ym(i) ,no(i) ,sg(i) , id(i) ,ary(i) ,anlyst(i) , 

&  gb(i) ,i*l ,ns) 

100  format  (i5,3f8.1,i5,2L5,i5) 
c 

c...  Main  call 

c 

cal 1  dynkin(ym,no ,ns , id , sg , fm, se , f 1 , ss , c , 

&  znop, zlrg, zv lrg, ary, anlyst, ymc, idc, pthr, pobs, gb,gbc, ok, 

&  pkin,naa,nad,nsa,nsd) 
c 

vrite(2,200)  ok 
write(2,z50)  fm,fl,pobs 

250  forma t('  fm  -  ',rl0.2,'  fl  -  ',  fl0.2/  pobs  -  ',fl0.2) 

200  format  ('  ok  •  ',L5) 

write(2,l00)  (i,ym(i) ,no(i) ,sg(i) ,id(i) ,ary(i) , anlyst (i) , 
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c 

c  •  • 
c  •  • 
c 

c .  . 
c. . 
c. . 
c. . 
c 

c. . 
c  •  » 

C  .  . 
c .  . 
c. . 
c 

c .  . 
c. . 
c. . 
c. . 
c. . 

C.  . 

c. . 
c. . 
c. . 
c. . 
c  •  • 
c  •  • 
c  •  • 
c. . 
c  .  • 
c 

c  •  « 
c  •  • 
c  •  • 
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c. . 
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C  •  • 
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c. . 
c  •  • 
c  •  • 
c 

C.  . 

c. . 
c. . 
c. . 
c. . 
c. . 
c. . 
c 


subroutine  dynkin(ym,no ,ns , id,sg,fm,se,fl,ss,c, 

&  znop, zlrg,zvlrg ,ary .anlyst , ymc, idc.pthr.pobs ,gb,gbc,ok, 

&  pkin,naa,nad,nsa,nsd) 

Subroutine  to  examine  amplitudes  of  an  event  and  determine 
if  they  are  anomalous. 

id(.)"l  if  observed  and  amplitude  measured 
id(.)“2  if  observed  but  amplitude  not  measured  so  signal 
greater  than  noise 

id(.)“3  if  not  observed  so  signal  is  less  than  noise 

idc(.)  same  as  id  except  with  changed  values  as  appropriate 
on  return 

gb(.)"l  if  station  is  good  and  should  not  be 
thrown  out 

gb(.)"0  if  station  is  bad  and  should  be  thrown  out 

gbc(.)  same  as  gb  except  with  changed  values  as  appropriate 
on  return 

ym  input  array  of  ns  magnitudes  with  distance  and  station 

corrections  applied 

ymc  returned  array  of  ns  magnitudes 

no  array  of  station  noise  levelB  with  distance  and 

station  corrections  applied. 

sg  noise  standard  error  (typically  0.2,  but  smaller 

if  noise  measured  especially  for  this  event) 
se  standard  error  of  network  magnitude  estimate 

fl  log-likelihood  of  fm 

as  standard  error  of  signal  amplitudes  (typically  0.4) 

might  want  to  make  this  an  array  if  signal  variance 
were  known  to  be  small  at  some  station 
c  S/N  required  for  detection 

pthr  probability  threshold  required  to  pass  kinematic 
test,  typically  .95. 

pobs  returned  kinemat^.  probability  of  the  event 
ok  ok». true.  ,  nothing  wrong  with  event  as  returned, 

if  ok  “  .false,  the  event  as  returned  fails  the 
kinematic  probability  test  and  should  be  released 
ary  ary  *  .true,  station  is  an  array  station 
anlyst  anlyst  “  .true,  detection  verified  by  analyst 
These  two  variables  are  used  to  define  four 
"types"  of  detections,  and  thus  the  "type"  event 
aa,ad  number  of  sites  with  array  and  analyst,  array 
and  automatic  detector. 

sa,sd  same  as  for  arrays  but  for  single  site 

snop  if  a  station  is  inoperative  it  will  incorrectly 

be  "reported"  as  not  detecting  because  the  signal 
is  below  the  noise.  If  this  station  is  left  out 
then  the  event  magnitude  will  increase,  i.e. 
z  will  be  negative,  znop  is  the  threshold  to 
detect  these  events;  we  should  set  its'  threshold 
fairly  high,  say  znop  "-.7 

C-l  7 
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c  •  •  • 
c  •  •  • 
c  •  •  • 

c 


c 


10 

c 

20 

c 

c  •  «  . 
c 


c 

c. .  . 
c 


c 

c. .  . 
c  •  •  • 
c 


zlrg  a  large  amplitude  generates  a  positive  z,  try 
zlrg  -  I 

zvlrg  z  for  throwing  out  wild  values,  try  zvlrg  “1-5 
Parameter  (DIMEN-50) 

Real  ym(ns ) ,ymc (ns ) ,no(ns ) ,sg(ns ) , fm, se , f 1 , ss , c , z(DIMEN) ,ss 

Integer  naa,nad,nsa,nsd 

Real  pkin ( 0 : na a , 0 : nad , 0 : n sa , 0 : n ad ) 

Real  zma, zmin, znop, zlrg , zvlrg ,pthr,pobs 
Integer  id (ns) ,idc(ns ) ,gb(ns ) ,gbc(ns ) ,ns , i 
Integer  izma, izmin.aa ,ad , sa , sd , izf 
Logical  ok,ary(ns),anlyst(ns) 

ok  “  .false* 
do  10  i  »  l,ns 
ymc(i)  ■  ym(i) 
idc(i)  ■  id(i) 
gbc(i)  *=  gb(i) 
continue 

call  type(ary>anlyst,idc,gbcln8,aa,ad,sa,sd) 

call  detz(ymc,no,ns,idc,sg,gbc,fm,8e,fl,ss,c,z) 
call  maxv(z,idc,gbc1ns, zma, izma) 
call  minv(z,idc,gbc,n8,zmin, izmin) 

first  throw  out  wild  ones 


if(abs(zma).gt. zvlrg. or. ab« (gmin). gt. zvlrg)  then 
if (abs(zma) .gt.abs(zmin))then 
izf- izma 


else 


else 

izf-izmin 

endif 


then  check  others 

it  ( zwin.gL . zuop. and . zma . iL . zlrg)  Lheu 
ok  ■  .true, 
return 

endif 

if  ( zmin. le. znop)  then 
izf*izmin 

else 

consider  large  amps  only  if  no  anomalous 
small  amplitudes  remain 

izf-izma 

endif 

endif 

if ((idc(izf ) .eq.3) .or . ( idc( izf ) .eq.2) )  then 
gbc(izf)  ■  0 
idc(izf)  *  3 

endif 

if(idc(izf ) .eq.l)  then 
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idc(izf)  =  2 
ymc(izf)  ■  0.0 

above  formula  assumes  noise  levels  are  in  magnitude 
units 

endif 

call  type (ary ,anlyst, idc,gb,ns ,aa,ad,sa,sd) 
if (aa. le.naa.and .ad. le .nad.and . sa. le.nsa.and . sd . le .nsd)then 
pobs  -  pkin(aa,ad,sa,sd) 

else 

pobs  *  .999 

endif 

if (pobs.gt.pthr)  then 
go  to  20 

else 


return 


returns  with  ok  =  .false 


endif 

end 


•  .*•  /•  .n  '  '.  v  , 
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subroutine  detz ( ym, no, ns, id , sg ,gb , fm, se , £1,88,0,2) 

Calculate  z(.),  an  array  of  normal  variables  showing 
influence  of  ym  measurement  on  fm.  After  Bob  Shumway 
June  1982 

id(.)=l  if  observed  and  amplitude  measured 
id(.)*2  if  observed  but  amplitude  not  measured  so  signal 
greater  than  noise 

id(.)*3  if  not  observed  so  signal  is  less  than  noise 
gb(.)=l  if  station  is  good  and  should  not  be  thrown  out 
gb(.)~0  if  station  is  bad  and  should  not  be  thrown  out 

ym  array  of  ns  magnitudes  with  distance  and  station 

corrections  applied 

sg  noise  standard  error  (typically  0.2,  but  smaller 

if  noise  measured  especially  for  this  event) 

"xxi"  variables  internal  to  amp  so  that  variables  don't 
scrambled 

se  standard  error  of  network  magnitude  estimate 

fl  log- likelihood  of  fm 

88  standard  error  of  signal  amplitudes  (typically  0.4) 

might  want  to  make  this  an  array  if  signal  variance 
were  known  to  be  small  at  some  station 
c  S/N  required  for  detection 

Parameter  (DIHEN  m  5 0) 

Real  ym(ns ) ,no (ns ) ,ag (ns ) ,se , f 1 , 88 , fm, z(ns ) 

Integer  id (ns) ,gb(na) ,ns 

Real  ymi(DIMEN) ,noi (DIMEN), sgi (DIMEN ) ,fmr,ser, fir ,c 
Integer  idi(DIMEN) ,ic,i,k, j,L,gbi (DIMEN) 

ic  ■  1 

call  amp(ym,no,n8,id,8g,gb,ic,fm,se,fl,ss,c) 
do  20  i  ■  l,ns 
k  -  1 

do  15  j  ■  l,ns 
if(j.eq.i)  go  to  13 
ymi(k)  -  ym(j) 
noi(k)  ■  no(j) 
idi(k)  ■  id(j) 
sgi(k)  -  sg(j) 
gbi(k)  -  gb(j) 
k  ■  k  ♦  1 
continue 

write(2,i01) (i,t,ymi(L) ,noi(L) , idi(L) ,L-1 ,ns-l) 
format(2i5,2fl0.3 ,i5) 

call  amp(ymi,noi,ns-l ,idi ,sgi ,gbi , ic,fmr,ser,f lr ,ss,c) 

z(i)"(fm-fmr )/ se 

continue 

return 
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Subroutine  amp(ym,no,n8,id,6g,gb,ic,fm,se,fL,6B,c) 

Calculate  log  likelihood  of  amplitude  configuration; 
or 

Maximum  likelihood  estimator  for  magnitude  using  Newton- 
Saphson  method 

by  Bob  Shumway,  June  1982 

ic  >0  if  likelihood  only  at  fm  needed 

ic  “1  if  Newton-Raphson  estimator  for  fm  needed 

id(.)*l  if  observed  and  amplitude  measured 
id(.)*2  if  observed  but  amplitude  not  measured  so  signal 
greater  than  noise 

id(.)*3  if  not  observed  so  signal  is  less  than  noise 

gb(.)*l  if  station  is  good  and  should  not  be  thrown  out 
gb (.)*0  if  station  is  bad  and  should  not  be  thrown  out 


array  of  magnitudes  with  distance  and  station 
corrections  applied 

array  of  noise  magnitudes  with  distance  and  station 
corrections  applied;  substituted  for  ym  when 
id.ne.l 

noise  standard  error  (typically  0.2,  but  smaller 
if  noise  measured  especially  for  this  event) 
standard  error  of  network  magnitude  estimate 
log-likelihood  of  fm 

standard  error  of  signal  amplitudes  (typically  0.4) 
might  want  to  make  this  an  array  if  signal  variance 
were  known  to  be  small  at  some  station 
S/M  required  for  detection 

magnitude  tolerance  for  terminating  iterations 


m 


toler 


Parameter  (toler  ■  .002) 

Integer  ns 

Real  ym(ns) ,no(ns ) ,sg(ns ) ,se, fL,ss, fm,c 

Integer  id(ns),gb(ns),ic 

Real  fmb,dl ,d2,z,ph,r ,w,88s, lc 

Integer  nl,it, j,k,nit 

Real  probf,p 

nit *20 
888*88*86 
lc  «  loglO(c) 
if(ic.eq.O)  go  to  15 
fmb^). 
nl-0 

do  10  j*l,ns 

if ((id( j) .gt.l) .or.(gb( j) .eq.0) )  go  to  10 
nl*nl+l 
fmb*f mb+ym( j ) 
continue 

£ob«fmb/ (Real (nl ) ) 
fm-fmb 
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go  to  17 
15  nit=l 
17  continue 

do  40  it  =  l,nit 
f  L=0 . 
dl-0. 
d2=0 . 

do  35  j=i,ns 
if (gb( j ) .eq .1 )  then 
k=id( j ) 

go  to  (20,25,30)  ,k 

else 

go  to  35 

end  if 

20  fL=f L-log(s8s)-( ( yml j )-fm)**2) / (2 .*ss  s) 

dl“dl+(fmb-fm)/ (sss) 
d2-d2-l ./ (sss) 
go  to  35 

25  w“8qrt(sss+8g( j)*sg(  j)) 
zm ( f  m-no ( j ) - lc ) / w 
ph“probf (z) 
r**p(z)/ph 
fL-fL+log(ph) 
dl“dl+r/w 

d2«d2-(z*r+r*r)/  (w*w) 
go  to  35 

30  w“sqrt(ssa+8g( j)*sg( j) ) 
z-(no( j)+lc-fm)/w 
ph-probf (z) 
r“p(z)/ph 
fL-fL+log(ph) 
dl*dl-r/w 

d2»d2+(-z*r-r*r )/ (w*v) 

35  continue 
c 

se“sqrt(-l ./d2) 

if ( it . eq .nit .or .abs (dl/d2) . It . toler )  go  to  40 
fm*fnrdl/d2 

c  write(2,l00j  it,tm 

100  format("  it,fm  , i5  , f 10 .3 ) 

40  continue 
return 
end 
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4.0 

3.y 

0.2 

1 .true . 

.true . 

3.6 

3.y 

0.2 

1 .true . 

.true 

4.4 

3.9 

0.2 

1  .false . 

.true 

2.9 

3.9 

0.2 

1  .false. 

.true 

0.0 

5.2 

0.2 

2. false. 

.true 

0.0 

2.0 

0.2 

3. false. 

.true 

0.0 

4.2 

0.2 

3  .false . 

.true 

0.0 

4.2 

0.2 

3  .false. 

.true 

0.0 

3.9 

0.2 

3. false. 

.true 

0.0 

4.5 

0.2 

3  .false. 

.true 

0.0 

5.0 

0.2 

3  .false. 

•  true 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 
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sta 

mag 

noise 

noi  var 

id 

ary  anlyst  gb 

1 

4.0 

3.9 

0.2 

1 

t 

t 

1 

2 

3.6 

3.9 

0.2 

1 

t 

t 

1 

3 

4.4 

3.9 

0.2 

1 

f 

t 

1 

4 

2.9 

3.9 

0.2 

1 

f 

t 

1  — 

5 

0. 

5.2 

0.2 

2 

f 

t 

1  — 

© 

6 

0. 

2.0 

0.2 

3 

f 

t 

1  — 

7 

0. 

4.2 

0.2 

3 

f 

t 

1 

8 

0. 

4.2 

0.2 

3 

f 

t 

1 

9 

0. 

3.9 

0.2 

3 

£ 

t 

1 

10 

0. 

4.5 

0.2 

3 

f 

t 

1 

11 

0. 

5.0 

0.2 

3 

f 

t 

1 

ok 

t 

fm  ■ 

3.88 

fl  - 

2.34  pobs  » 

0.95 

1 

4.0 

3.9 

0.2 

1 

t 

t 

1 

2 

3  .6 

3.9 

0.2 

1 

t 

t 

1 

3 

4.4 

3.9 

0.2 

1 

f 

t 

1 

4 

2.9 

3.9 

0.2 

1 

f 

t 

1 

5 

0. 

5.2 

0.2 

2 

f 

t 

1 

6 

0. 

2.0 

0.2 

3 

f 

t 

1 

7 

0. 

4.2 

0.2 

3 

f 

t 

1 

8 

0. 

4.2 

0.2 

3 

f 

t 

1 

9 

0. 

3.9 

0.2 

3 
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t 

1 

10 

0. 

4.5 

0.2 

3 

f 

t 

1 

11 

0. 

5.0 

0.2 

3 

f 

t 

1 

sta 

mag 

idc 

gb 

1 

4.00 

1 

1 

2 

3.60 

1 

1 

3 

4.40 

1 

1 

4 

0. 

2 

1 

5 

0. 

3 

0 

6 

0. 

3 

0 

7 

0. 

3 

1 

8 

0. 

3 

1 

9 

0. 

3 

1 

10 

0. 

3 

1 

11 

0. 

3 

1 

Detected,  noise  reasonable  for  detection  but  signal  measured  too 
low,  so  retain  detection  but  discard  amplitude  and  say  9 >  N,  so  tdc 
goes  1  — *•  2 

Detected,  but  noise  too  high,  so  discard  detection  -  e.g.  a  shadow 
zone  "detection".  So  gb  goes  1 — *-0 

©  Mot  detected,  but  noise  very  low  -  so  discard  station,  it  must  be 
down  so  gb  goes  l  — *•  0 
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subroutine  event (ym, no, ns , id, sg ,ary .anlyst ,gb) 

Integer  DIMEN 

Parameter  (DIMEN  *  501 

Integer  id(DIMEN) , ns, i,gb( DIMEN) 

Real  ym(DIMEN) ,no(DIMEN) ,sg(DIMEN) 

Logical  ary (DIMEN) , anlyst (DIMEN) 

Sample  data  from  8  September  82  memo  by  Blandford 
and  Shumway.  Would  typically  be  replaced  by  access 
package  to  files. 

Do  10  i-1, DIMEN 
ym(i)*0. 
sg(i)«0. 
id(i)*0 

ary (i)«. false, 
anlyst (i)s. false. 
gb(i)*0 
)  continue 
nsBll 

open(l, status* 'old' ,formF 'formatted' ,file*'event. in' ) 
rewind  1 

read  (1,1)  (ym(i) ,no(i) ,sg(i) ,id(i) , ary (i) .anlyst (i) ,gb(i) , 
&  i*l,ns) 

format (3f5 .1 ,i5 ,2L10,i5) 
return 
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c 

c «  •  • 
c.  •  . 


c  •  •  • 


c 


10 


subroutine  maxv(z,id,gb,ns,zma,izma) 

Real  z(ns),zma 

Integer  ns, izma , i , id(ns ) ,gb(ns ) 

Finds  max  value  of  z  statistics  for  good  stations 
Initial  value  of  zma  is  only  -10  because  are 
working  with  standard  deviation  units 


zma *-10 . 
do  10  i  »  l,ns 

if  (z(i)  .gt.zma.and.gb(i)  .eq.Dthen 
zma  “  z(i) 


end  if 
continue 
return 
end 


izma 
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subroutine  minv(z, id,gb ,ns , zmin, izmin ) 
c 

c...  gets  maximum  of  z  statistics  for  good  stations 

c...  zmin  initially  so  small  since  working  in  standard  deviations 

c 

Real  z(ns),zmin 

Integer  ns,izmin,i(id(ns),gb(ns) 
zmin  ■  10. 
do  10  i  ■  l,ns 

if(z(i).lt.zmin.and.gb(i)  .eq.Dthen 
zmin-z(i) 
izmin  ■  i 

endif 

10  continue 
return 
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subrout ine  mkpkin ( pkin , nnaa, nnad ,nn sa , nn  sd ) 
c 


c. 

•  • 

Reads 

in  the 

kinematic 

probability  matrix. 

c. 

e  • 

This 

assumes 

that  Nad 

and 

c. 

•  • 

Nsd  are  0.  This  is  an 

example  only; 

the  real 

matrix 

should 

c. 

c . 

•  • 

be  derived  from  analysis  of  true 

and 

false  events. 

c. 

•  • 

nsa 

c. 

•  e 

0 

1 

2 

3 

4 

5 

6 

7 

c 

naa 

c 

5 

.999 

.999 

.999 

etc. 

c 

4 

.99 

.99 

.999 

etc . 

c 

3 

.97 

.97 

.97 

.99 

.999 

etc. 

c 

2 

.9 

.95 

.95 

.99 

.999 

etc 

c 

1 

.0 

.5 

.8 

.95 

.97 

.99 

.99 

.999 

c 

0 

.0 

.0 

.0 

.2 

.8 

.9 

.95 

.97 

c 

Integer  naa,nad,nsa,nsd,nnaa,nnad,nnsaannsd 
Real  pkin(0:nnaa,0:nnad,0;nnsa,0:nnsd) 
open(3,  status*' old'  .form* 'formatted' ,  file*  'pkin.  in') 
rewind  3 

read(3 ,100) ( ( ((pkin(naa,nad,nsa,nsd ) ,naa*0 ,nnaa) ,nad*0 ,nnad ) , 
&  nsa*0,nnsa) ,nad*0,nnsd) 

100  foraat(6f5.3) 
return 
end 
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function  p(x) 

Real  p,x 

p=exp(-x*x/2.)/2. 50662828 

return 

end 


029 
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function  probf(x) 
double  precision  q,t 
Real  x,probf 
t*abs(x) 

q-(( ((t*. 0006341 9-. 00010754)*t+. 01057 706)*t+.04833145)*t+ 
&  .10882473)*t+l .09050773 

q-1./ (q*q) 

q-q*q 

probf«q*q 

if(x.le.0.)  go  to  2 

1  probf*l .-probf 

2  return 
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subroutine  type(ary ,anlyst, id,gb,ns,aa,ad,sa,sd) 

Integer  ns, id (ns ) ,gb(ns  )  >aa , ad ,  sa , sd  ,  i 

Logical  ary(ns) ,anlyst(ns) 

aa=0 

ad“0 

sa“0 

sd*0 

do  10  i  ■  l(ns 

if(gb(i).eq.0.or.id(i).gt.2)  go  to  10 
if(ary(i)«and.anlyst(i))  aa  *  aa  +  1 

if (ary(i) .and . .not.anlyst (i) )  ad  =  ad  +  1 

if ( .not.ary(i) .and.anlyst(i))  sa  ■  sa  +  1 

if ( .not.ary(i) .and. .not.anlyst(i))  sd  ■  sd  +  1 

10  continue 
return 


MEMORANDUM 


m 


FROM: 


DATE: 


J.  Goncz 


R.  Blandford  'i-  * 


SUBJECT:  A  Better  Approach  to  Handling  Amplitude  Data  in  AA 


19  April  1983 


1.  Compute  maximum  likelihood  m,  using  the  observing  stations  plus  a 
fixed  set  of  about  20  reliable  stations.  An  is  calculated  with 
each  station  omitted  in  turn.  Assume  an  a-priori  a  *  0.35  m^ 
throughout . 

2.  Analyze  each  of  the  above  stations  in  screen  as  we  do  now,  get  the 
total  number  of  points,  then: 

i.  If  the  station  did  not  detect  but 


p(det)  >  .95,  -1  point 
p(det)  >  .99,  -2  points 

Comments:  If  p(det)  >  .95  at  many  stations  most  stations  will 

have  detected  if  the  event  is  big  and  loss  of  points  for  the 
occasional  random  non-detection  will  not  matter.  If  the  event 
is  weak,  only  a  few  close  stations  will  have  p(det)  >  0.95  so 
the  occasional  random  non-detection  is  again  unlikely  to  lose 
points . 

If  the  station  did  detect  and  if,  based  on  noise  and  m^  alone 

p(det)  <  .01, -(all  points  due  to  station) 

Comments:  If  the  event  is  big  there  will  be  so  many  points  that 

occasionally  losing  one  will  not  be  serious.  If  it  is  a  small 
event  and  just  trapped  (correctly)  the  one  distant  station,  then 
the  probability  of  that  happening  for  any  particular  one  out  of 
the  20  distant  stations  may  be  smaller  than  .05  so  for  that 
reason  we  set  this  threshold  at  .01  so  that  only  about  1/5  of 
these  events  will  be  incorrectly  rejected.  It  is  important  to 
reject  all  the  points  due  to  the  station,  i.e.  those  due  to 
azimuth  and  slowness  and  those  due  to  local  flags  and  associated 
S  phases. 

If  the  station  did  detect,  and  reports  an  amplitude  Aq,  and  if, 
based  on  noise  and  m^  alone, 

p(det)  >  .01  and 

p(amp  >  Aq)  <  .01  or 

p(amp  <  A  )  <.01  then 

o 
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discard  amplitude,  A  (are  left  with  estimated  noise)  and  return 

to(l)  0 

Comments:  That  is,  compute  all  the  negative  points  and  m, 

values  over  because  the  bad  amplitude  probably  distorted  all  the 
m^  values.  Note  that  this  procedure  does  not  "fiddle"  with  an 
event.  It  either  results  in  complete  acceptance  or  discards  the 
event  so  that  the  kinematic  processing  has  a  new  chance  to  do 
better. 


1.  Steps  Taken  in  the  Generation  of  Short-Period  Arrivals 


(1)  For  each  event-station  pair,  a  true  event  distance  was  computed 

(2)  For  each  phase  valid  at  this  distance,  B-factor  tables  were  used  to  compute 
a  station  value  for  the  amplitude /period  ratio  of  the  arrival  For  P,  distance 
<100  degrees,  the  B-factors  given  by  Veith  and  Clawson  (1969)  were  used  . 
for  all  other  phases,  the  B-factors  given  in  the  tables  at  the  end  of  this 
appendix  were  applied.  The  mb  bias  values  listed  in  table  1  of  USVGSF/7 
were  included  for  P  only.  Scattering  in  the  earth  was  modelled  by  a  Gaus¬ 
sian  term  with  standard  deviation  0.3  mb  units.  The  period  T  was  then 
selected  randomly  in  a  small  interval  about  1  second  and  the  amplitude  A 
computed  in  nanometers. 

(3)  For  events  at  distances  D  less  them  6  degrees,  Pg  amplitudes  A  were  com¬ 
puted  using 

log(A)  =  mb  -  1.6  -  3  log(D) 
and  S  amplitudes  were  taken  to  be  twice  those  of  P. 

(4)  The  amplitude  A  of  each  phase  was  then  compared  to  the  station  noise 
values  given  in  CCD/558.  If  the  signal  to  noise  ratio  was  less  than  0,  ihc 
arrival  was  not  accepted;  if  between  1  and  2  ,  it  was  accepted  with  a  random 
probability  of  50%  ;  if  between  2  and  3,  it  was  accepted  with  a  random  pro¬ 
bability  of  75%;  if  greater  than  3,  the  arrival  was  accepted 

It  should  again  be  noted  here,  as  in  Chapter  2,  that  the  noise  levels  given  in 
CCD/558  were  doubled  in  order  to  keep  the  number  of  arrivals  within  rea¬ 
sonable  bounds.  A  Gaussian  term  with  0.2  mb  unit  variance  was  added  to 
the  noise  before  computing  the  signal  to  noise  ratio 

(5)  Station  downtime  was  modelled  by  allowing  the  user  or  the  program  to 
specify  two  down  intervals  per  day  for  each  station. 

(8)  The  phase  type  of  arrivals  reported  was  always  given  as  P  or  blank 
(unnamed)  according  to  (13)  below,  except  for  depth  phases  (see  (11) 
below)  and  Pg  and  S  phases,  which  were  correctly  named 

(7)  The  standard  deviation  of  travel  times  was  set  at  1  second  for  all  phases  and 
distances. 

(8)  Arrival  times  of  the  onset  were  then  computed  from  the  true  event  time 

(9)  To  simulate  operator  error,  in  1%  of  all  arrivals  (selected  randomly)  an 
error  of  1  was  included  in  the  minute  designator  (randomly  positive  or 
negative).  A  10-second  error  was  randomly  introduced  in  3%  of  all  cases 
except  for  array  stations. 

(10)  For  array  observations,  an  error  in  vector  slowness  with  standard  deviation 
1.5  seconds /degree  was  added  to  each  component  of  the  slowness  vector 
calculated  from  the  true  azimuth  and  tabulated  slowness  for  each  phase 

(11)  Depth  phases  are  named  according  to  the  following  rules 

a.  If  P  is  not  detected,  pP  or  sP  are  not  named  (blank). 

b.  If  P  and  sP  are  detected,  but  pP  is  not,  and  (P-sP)time  <i  min  ,  label  sP  as 

pP. 

c.  If  P  and  PcP  are  detected  but  pP  and  sP  are  not,  and  (P-PcP)timc  <  1  mm., 

label  PcP  as  pP 


(12)  After  all  arrivals  for  a  given  station  and  tune  interval  have  been  generated, 
the  following  process  is  applied  to  screen  out  arrivals  that  probably  could 
not  have  been  determined  because  of  interference  from  larger  earlier 
arrivals 

a.  As  a  definition,  there  is  a  signal  S  with  amplitude  As,  possibly  hidden  by  a 
previous  signed  of  amplitude  A. 

b.  !f  there  is  not  a  signed  in  the  60  seconds  before  S,  let  S  exist. 

c.  If  there  is  a  signal  in  the  previous  60  seconds,  let  S  exist  if  As  >0.1A,  other¬ 
wise  reject  S. 

d.  If  there  is  a  signal  in  the  previous  10  seconds,  let  S  exist  if  As  >0.5A,  other¬ 
wise  reject  S. 

e.  Even  if  S  is  rejected,  it  can  still  hide  arrivals  coming  after  it  according  to  c. 
and  d.. 

(13)  Except  for  Pg  and  S  (see  (6)  above  )  and  depth  phases  (see  (ll)  above), 

phases  are  named  as  follows  : 

a.  If  P,  named  P  unless  there  is  another  arrival  less  than  1  minute  beforehand, 
in  which  case  it  is  unnamed. 

b.  If  anything  other  than  P.  named  P  if  there  have  beenJarrivals  in  the  preced¬ 
ing  3  minutes  ;  otherwise  unnamed. 
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B-FACTORS  for  SECONDARY  PHASES 


Phase 


Dist. 

B 

Phase 

Dist. 

B 

h<70 

B 

h>70 

10 

3.0 

SKPdf 

105 

5.4 

4.6 

11 

3.1 

110 

5.4 

4.6 

12 

3.2 

113 

5.2 

4.3 

13 

3.3 

117 

5.0 

4.1 

14 

3.4 

120 

4.8 

3.9 

15 

3.5 

125 

4.7 

3.8 

20 

3.7 

130 

4.6 

3.7 

25 

3.9 

140 

4.4 

3.5 

30 

4.0 

180 

4.4 

3.5 

35 

3.7 

SKPab 

all 

3.9 

3.7 

40 

3.6 

ScP 

all 

4.4 

3.9 

45 

3.6 

pP 

all 

B(P) 

B(P) 

50 

3.7 

-0.3 

-0.5 

55 

3.8 

sP 

all 

B(P) 

B(P) 

60 

3.9 

-0.3 

-0.5 

HP 


PcP 

PKKPab 

PKKPbc 

PKKPdf 


PKPPKP 


65 

70 

75 

80 

>80 

0-80 

all 

all 

20 

30 

40 

50 

60 

130 

140 

150 

160 

170 

180 

30 

40 

50 

55 

95 

100 

110 

120 

130 

140 


4.0 

4.1 

4.1 

4.2 

4.2 
4.00 
4.60 
4.60 
5.5 
5.4 

5.2 
5.0 
4.7 

4.7 

4.8 
5.0 

5.2 

5.4 

5.5 
5.4 

5.3 
5.0 

4.6 
4.6 
5.0 

5.3 

5.4 
5.4 
5.4 
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